Asymptotic decay of pair correlations in a Yukawa fluid 



in 
o 
o 

(N 
o 

Q 



p. Hopkins, A.J. ArcherQ and R. Evans 
H.H. Wills Physics Laboratory, University of Bristol, Bristol BS8 ITL, United Kingdom 

(Dated: February 2, 2008) 

We analyse the r — > cxd asymptotic decay of the total correlation function, h{r), for a fluid 
composed of particles interacting via a (point) Yukawa pair potential. Such a potential provides 
a simple model for dusty plasmas. The asymptotic decay is determined by the poles of the liquid 
structure factor in the complex plane. We use the hypernetted-chain closure to the Ornstein-Zernike 
equation to determine the line in the phase diagram, well-removed from the freezing transition line, 
where crossover occurs in the ultimate decay of h{r), from monotonic to damped oscillatory. We 
show: i) crossover takes place via the same mechanism (coalescence of imaginary poles) as in the 
classical one-component plasma and in other models of Coulomb fluids and ii) leading-order pole 
contributions provide an accurate description of h{r) at intermediate distances r as well as at long 
range. 
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The Yukawa, or screened Coulomb potential, is often 
used as a model for fluids composed of charged parti- 
cles immersed in a uniform neutralising background. We 
define the pair potential as 



to be nil 



eexp(— Ar) 



(1) 



where r is the distance between the centres of a pair of 
particles, A is the inverse decay length (screening parame- 
ter) and e > characterises the strength of the potential. 
A fluid composed of particles interacting via this poten- 
tial exhibits behaviour similar to that of the one compo- 
nent classical plasma (OCP) [3,13, in which (j>{r) cx l/r. 
The point Yukawa pair potential also corresponds to a 
limiting case of the well established Derjaguin, Landau, 
Verwey and Overbeek (DLVO) model for aqueous sus- 
pensions of charged colloidal particles H, H, |^ , where the 
hard-core part of the colloid-colloid effective pair poten- 
tial is neglected. 

The Yukawa pair potential is also widely employed in 
theoretical studies of the so-called dusty plasmas. These 
are multicomponent plasmas consisting of charged (dust) 
particles, electrons and ions, as well as neutral atoms or 
molecules, which are found in a variety of environments, 
from the interstellar medium to plasma etching processes. 
Depending on their size, the dust grains can attain a 
large negative charge of 1000-lOOOOe for particles of size 
1-10/iTO 1^; the charge is generally negative and is de- 
termined by the balance of the absorbed electron and 
ion fluxes. Since the dust component of the plasma can 
be videoed and tracked directly, dusty plasmas provide a 
valuable system for studying both equilibrium phenom- 
ena and collective processes such as transport in a fluid 
0, Q ■ In modelling the dusty plasma, the effective po- 
tential between two dust particles of charge Q is taken 



(2) 



where is the Debye length of the background plasma. 
The thermodynamics of the system can then be charac- 
terised by the dimensionless parameters 



k]ja and F 



AneoaksT 



(3) 
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where a = (3/(47rp))^/^ is the Wigner-Scitz radius, i.e. 
the mean interparticle distance, which is determined by 
the average fluid number density p. T is the coupling 
(plasma) parameter Note that unlike the OCP, whose 
properties depend solely on F, two parameters are re- 
quired in this case. 

Because of its relevance for colloidal fluids and the 
growing importance of studies of dusty plasmas, not to 
mention the appealing mathematical properties of the 
Yukawa potential ^3 > much is established for this model 
fluid. The phase diagram has been determined in a num- 
ber of simulation studies in El III 111 III III At small 
values of F the liquid state equilibrium structure is well- 
approximated by the hypernetted-chain (HNC) closure 
to the Ornstein-Zernike (OZ) equation ll^ ITsl |T9|. as 
one would expect from studies of the OCP jjl- For larger 
values of F, a modified HNC based upon a re-scaling of 
the bridge function for the OCP |22| yields results in- 
distinguishable from Monte-Carlo (MC) simulation data 
. At small values of k the Yukawa fluid freezes into a 
bcc crystal upon increasing F, as in the OCP. For suffi- 
ciently large k increasing F can lead to freezing directly 
into a fee crystal. A portion of the phase diagram from 
Ref. nil is shown in Fig.[T] llj . 

In the present work we analyse the asymptotic decay 
of the pair correlations in the uniform fluid and we find 
that at couplings F below freezing, there is a crossover 
in the form of the asymptotic decay, r ^ cxo, of h{r), 
the total correlation function, from monotonic to damped 
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FIG. 1: Crossover line (diamonds joined by a solid line) 
separating the region of the (k, F) plane where the asymp- 
totic decay of h{r) is damped oscillatory from that where it 
is monotonic. In the OCP, corresponding to k = 0, crossover 
occurs at Tk ~ 1.21 f22| . We also display the fluid-solid (-f) 
and solid-solid (x) phase boundaries given in Ref. p^ . In the 
inset we display the crossover line in the (p, T) plane. Note 
that the freezing transition present at low T* = ksT/e is not 
visible on this scale. 



oscillatory, that is similar to the crossover found near 
F/f = 1.12 in the HNC treatment of the OCP jl^. We 
map-out the crossover line in the (k, F) phase diagram - 
see Fig. ^ 

Our starting point is the OZ equation which re- 
lates h{r) to c(r), the pair direct correlation function. In 
Fourier space the OZ equation can be written as: 



h{q) 



1 - pc{q) 



(4) 



Here f{q) denotes the 3-dimensional Fourier transform 
of the spherically symmetric function /(r). We choose to 
implement the HNC closure which sets the bridge func- 
tion B(r) — 0. For the present Yukawa fluid, particularly 
at the densities in the neighbourhood of the crossover 
line, the HNC accounts very well for the structure of 
the uniform fluid, yielding pair correlation functions al- 
most indistinguishable from simulation results. Thus we 
shall use the HNC results to determine the location of 
the crossover line. First in Fig. [3 a) we display some 
typical results for g{r) ~ 1 + h{r). These are for a 
flxed temperature T* = ksT/e — 1 and for the densities 
pX^^ = 0.1, 0.45, 2.5, 10 and 40. At this temperature, 
we shall find that the crossover occurs for pX~^ ~ 0.47. 
We also display some constant NVT MC results for the 
first three of these densities; for p\^^ ~ 0.1 and 0.45 we 
used 300 particles with 25000 (plus 10000 equihbration) 
trial moves per particle. The maximum displacement of 
the particles was chosen so that approximately 50% of at- 
tempted moves were accepted. For the MC simulations 



at p\~^ = 2.5 we used 2000 particles. We see as the 
density is increased the asymptotic decay switches from 
monotonic to damped oscillatory. This crossover repre- 
sents the evolution of the system from a weakly coupled 
state where the particles do not order strongly, to a state 
where they are more closely packed and the correlations 
become more hard-sphere like, although for all the densi- 
ties displayed, the amplitude of the (oscillatory) structure 
in g(r) remains quite small. 

The asymptotic decay of h{r) can be obtained [2^12^ 
from the OZ equation using the inverse Fourier transform 
of Eq. Q: 



rh{r) 



1 



dq q exp(igr) 



1 - pc{q) ' 



(5) 



which can be transformed into a semi-circular contour in- 
tegral in the upper half of the complex plane. Its value is 
determined by the poles of h{q) enclosed. These occur at 
qn — ±ai + ioif), where q^ is the solution to the equation 

1 - pc(g„) = 0. (6) 
As a result, h{r) can be obtained as the sum p^. l2^ 

"rHr) = 77- V'-RT! exp(zg„r), (7) 

n 

where q„ is the nth pole and i?„ is the residue of 
qc{q) / (l ~ pc(q)) at q„. Clearly the asymptotic behaviour 
of h{r) is determined by the pole(s) with the smallest 
imaginary part, oq. If this pole is purely imaginary, 
qn = ictQ, then rh(r) ~ ylexp(— ao^), for r 00, where 
yl is a (real) amplitude 23] . Alternatively, if the conju- 
gate pair qn = ±ai +iao has the smallest imaginary part 
then rh(r) ~ j4exp(— agr) cos(air — 9). The amplitudes 
A and A and the phase 9 can be calculated from the 
residues |23j . Whether a pure imaginary or complex pole 
dominates depends on the thermodynamic state point. 

In order to calculate the poles, we use the separation 
method introduced in Ref. Owing to the partic- 

ular form of the decay of the Yukawa pair potential, 
the asymptotic behaviour of c(r) must be treated sep- 
arately so as to ensure the convergence of the integrals 
which determine the poles The asymptotic decay, 

r — > 00, of the direct correlation function is given by 
c(r) ^ —(3(f){r), which for the Yukawa potential implies 
c(r) ~ — exp(— Ar)/(T*Ar). It is convenient to define a 
short-ranged direct correlation function c*'"(r) by sub- 
tracting the long-ranged Yukawa decay. The Fourier 
transform of c(r) is then 



47r 



1 



AT* (g2 + A2) ■ 



(8) 



Making this division, we follow Ref. |23| and calculate the 
poles by separating Eq. ^ into its real and imaginary 
parts and solving numerically using a Newton-Raphson 
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FIG. 2: Pair correlation functions for a reduced temperature 
T* = 1. a) g{r) calculated for a range of densities; the lines 
denote HNC results and the symbols MC results (see legend) . 
Crossover occurs at pX~^ — 0.47 for this temperature. For 
densities larger than this the decay of g{r) is damped oscil- 
latory, b) Comparison between the full HNC result for h{r) 
(solid line) and the contribution from the leading order sin- 
gle imaginary pole (dashed line) for pX'^ = 0.1. c) The fuU 
HNC result for h{r) (solid line) and the contributions from the 
two leading order purely imaginary poles (dotted and dashed 
lines) and their sum (dot-dashed line), for pX~^ — 0.45. d) 
The full HNC result for h{r) (solid line) and the contribu- 
tion from the leading order conjugate pair of complex poles 
(dotted line), giving damped oscillatory decay, for pX~^ = 40. 
The insets in c) and d) show ln[r-|/i(r-) |] versus Xr. 



FIG. 3: The leading order poles (smallest imaginary part, ao) 
along the isotherm T* = 1 for increasing pX~^. Only poles 
with positive ai are shown; complex poles occur in conjugate 
pairs. At low densities the leading order pole (x) is purely 
imaginary. As the density is increased a second purely imag- 
inary pole {+) descends and, at a density pA~^ ~ 0.47, the 
two purely imaginary poles coalesce to form a pair of complex 
poles (□). At this (Kirkwood) point the asymptotic decay of 
h{r) crosses over from monotonic to damped oscillatory. 



procedure. However, in general the integrals involved 
converge only for complex q such that ^[g] < 2ao, where 
ao is the imaginary part of the leading order pole, i.e. 
that with the smallest value of ao. In practice, the other 
poles generally lie outside this region of convergence and 
only the leading order pole can be determined. We also 
use this separation of c(r) to calculate the amplitude and 
phase of h{r) from the residues of the poles, assuming 
these to be simple p^ . 

Using our HNC results for c(r) we were able to cal- 
culate the contributions to h{r) from the leading order 
pole(s) for various points in the phase diagram. In Fig. [21 
b) we display the HNC result for h{r) at the state point 
pX^^ =0.1 and T* — 1, together with the contribution 
from the leading order (purely imaginary) pole. This 
state point lies on the monotonic side of the crossover 
line. In Fig.[5]c) we display the HNC result for h{r) at 
pX^^ = 0.45 and T* = 1 (near the crossover line, still 
on the monotonic side), together with the contributions 
from the two leading order (purely imaginary) poles. In 
Fig-Eld) we display the HNC result for h{r), for the state 
point pX~'^ = 40 and T* — 1, which exhibits damped 
oscillatory asymptotic decay, together with the contri- 
bution from the leading order conjugate pair of complex 
poles. Note that the contribution from this pair of poles 
approximates accurately the full h{r) for Xr > 0.2. 

In Fig. Owe display the leading order poles calculated 
along the isotherm T* — 1. At low densities the pole 
dominating the decay of h{r) is purely imaginary. How- 
ever, as the density is increased, this pole moves up the 
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imaginary axis, and at a density p\~^ — 0.47 it meets a 
second (descending) purely imaginary pole. This second 
pole could only be determined once it had descended into 
the region of convergence. These poles coalesce and form 
a conjugate pair of complex poles which then move away 
from the imaginary axis as the density is increased fur- 
ther. This coalescence of two purely imaginary poles to 
form a pair of complex poles as one moves along a path 
of increasing density in the phase diagram, results in a 
crossover in the asymptotic decay of h{r) from monotonic 
to damped oscillatory. The mechanism of poles coalesc- 
ing and moving off the imaginary axis is the same as 
that found in the OCP near Tk = 1.12 It is equiva- 
lent to that found in studies of charge correlations in the 
restricted primitive model (RPM) of binary ionic fluids 
|24| and in screened versions of the RPM 25] . Kirkwood 
[23 was the first to describe the mechanism, so the line 
in the phase diagram at which crossover occurs could be 
termed the Kirkwood line following earlier terminology 
1.22, 24, 25]. 

In Fig.^we display the location of the crossover (Kirk- 
wood) line in the phase diagram. This line (see inset) was 
determined by finding the density at which the imagi- 
nary poles coalesce for a series of isotherms, i.e. fixed T* . 
By converting from the variables {pX~^, T*) to (k, P) the 
crossover line shown in the main figure was obtained. 
We emphasise that this line lies far below the fluid-solid 
transition line, as might have been expected from the ob- 
servation that in the OCP freezing is known to occur for 
F ~ 172, whereas crossover from monotonic to oscilla- 
tory decay of correlations occurs near Tk ~ 1.12. For 



the values of k considered here, < k < 5, crossover 
occurs for F in the range 1 ^ F ^ 5, i.e. the coupling 
parameter is rather weak and we expect the HNC to per- 
form rather well; see also Fig. [21 a) for the comparison 
with MC. Of course, one could attempt to improve upon 
the HNC by incorporating an approximate bridge func- 
tion -B(r), along the lines of Refs. 0,0]. However, we 
do not expect significant changes in the location of the 
crossover line; the poles should not be sensitive to details 
of the particular closure [23] ■ One could also attempt to 
extract the poles and determine the crossover line using 
accurate simulation data for g{r) following the method 
applied in Ref. ;27] for a truncated Lennard- Jones fluid. 

In summary, we have shown that the form of the 
asymptotic decay of the total correlation function h{r) 
in a Yukawa fluid crosses over from monotonic at small 
coupling parameter F to damped oscillatory at larger 
values via the same mechanism as in the OCP. We 
find that leading-order asymptotics provide as accurate 
a description of pair correlations at intermediate range 
in the Yukawa fluid as they do for other model fluids 
mmil mill]. T^s observation might prove use- 
ful in further applications of the Yukawa model to dusty 
plasmas. 
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